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Abstract 

In this paper we present numerical simulations of the evolution of plan- 
ets or massive satellites captured in the 2/1 and 3/1 resonances, under 
the action of an anti-dissipative tidal force. The evolution of resonant 
trapped bodies show a richness of solutions: librations around station- 
ary symmetric solutions with aligned periapses (A-ro = 0) or anti-aligned 
periapses {Azu — tt), librations around stationary asymmetric solutions 
in which the periapses configuration is fixed, but with Auj taking values 
in a wide range of angles. Many of these solutions exist for large values 
of the eccentricities and, during the semimajor axes drift, the solutions 
show turnabouts from one configuration to another. The presented results 
are valid for other non-conservative forces leading to adiabatic convergent 
migration and capture into one of these resonances. 

Keywords: exoplanets, Galilean satellites, resonance, planetary mi- 
gration, tidal forces. 

1 Introduction 

The present paper considers the problem of a pair of planets (or planetary 
satellites) captured in the 2/1 and 3/1 mean-motion resonances, in which the 
inner mass is under the action of an anti-dissipative force that continuously 
increases its semimajor axis. Convergent migration of bodies and the action 
of non-conservative disturbing forces are subjects that have been explored by 
many authors (for a recent publication and bibliographical list, see Jancart 
2002). However, in almost all those studies, one of the bodies is assumed to 
have negligible mass and does not exert any gravitational action on the rest of 
the system. The scenario usually considered is thus the restricted three-body 
problem to which a non-conservative force is added. In this work, we are in- 
terested in systems where the two planets (or satellites) have finite masses. As 
examples of such systems, we may cite some of the newly discovered extraso- 
lar planetary systems (e.g. Gliese 876, HD82943, 55 Cue) and, in our Solar 
System, some satellite pairs such as lo-Europa or Europa-Ganymede. Even if 
we expect a phenomenology showing many of the characteristics seen in the re- 
stricted models, only actual calculations can tell us how much of the dynamics 
is equivalent and, if not, which are the main changes. 

The paradigm of a system of finite mass bodies captured into resonance is 
the system formed by the three inner Galilean satellites of Jupiter. If we denote 
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by Ai, A2 and A3, the mean longitude of lo. Europa and Ganymede, respectively, 
it was shown by Laplace that the motion of these bodies is such that the angle 

Ai — 3A2 + 2A3 

oscillates around tt with a small amplitude, ^ 10""^ rad (Lieske, 1998). Laplace 
also investigated the effect of a perturbation leading to an acceleration of the 
mean longitudes of one of the satellites. He showed that the mutual interaction 
of the three satellites would redistribute the effect among all the longitudes 
in such a way that the libration of the angle Ai — 3A2 + 2A3 around n would 
be conserved (see Tisserand 1896, Ferraz-Mello 1979). Nonetheless, the current 
tidal acceleration of the Galilean satellites is very small. Modern determinations 
of the rate acceleration/velocity, as deduced from old satellite eclipses and recent 
mutual events observations, are (2.88 ±0.10) x 10~^°yr~^ (Aksnes and Franklin 
2000) and (-0.074±0.087) x lO^^V"^ (Lieske 1987). Heat dissipation observed 
on the surfaces of lo and Europa indicates that tides raised by Jupiter on the 
satellites are among the main dissipative agents at work in this system (see 
Yoder, 1979). 

More recently, several resonant extrasolar planetary systems have been dis- 
covered with dynamical features suggesting that they have undergone a large- 
scale migration. In most cases, this orbital migration is believed to be due to an 
interaction of the planets with a residual planetesimals disk. In some extreme 
cases, tides have also been considered among the sources of orbital evolution; 
this seems to be the case of the putative planetary system around HD 83443 
(Wu and Goldreich, 2002) where the inner planet has the smallest semimajor 
axis among all currently known exoplanets. However, recent observations have 
not confirmed the existence of a pair of planets around that star. It is worth 
noting that many of the discovered exoplanets, including the inner planet of 
wAnd, orbit at distances smaller than 0.077 AU where tidal interaction with 
the central star should be important. All these planets are candidates to have 
an initial dynamical evolution similar to that shown by the Jupitcr-Io system. 

In this paper, we present a series of numerical simulations of the dynamical 
evolution of planetary satellites due to tidal effects, concentrating on the orbital 
variation after a capture in the 2/1 and 3/1 mean-motion resonances. Except 
for a change in the timescale of the process, these results should also be repre- 
sentative of the evolution of extrasolar planets moving very close to the central 
star. 

2 2-body tidal interaction 

The anti-dissipative force used in this paper corresponds to tides raised on the 
central mass mo (star or planet) by an orbiting body mi (planet or satellite). In 
addition, mi is interacting with a similar body m2 placed in an orbit external to 
that of nil. The model is planar and the spin axis of the central mass is assumed 
normal to the orbital plane. The two secondary bodies have finite masses but 
their sizes are neglected; they are thus considered as mass points. We will 
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further assume that the tidal interaction only affects the inner body and not 
the outer one. This is, by far, our most stringent hypothesis. Tidal forces are 
inversely proportional to the 7*^ power of the distance to the central body; at the 
beginning of our simulations, we always considered a\l < 0.6, making the tidal 
interaction of the outermost orbiting body lO'^ times smaller, but in simulations 
where the bodies escape an early capture and ai/a2 grows significantly, the 
results of a simulation taking into account both tidal interactions may be very 
different. 

Before presenting the physics of the used model, we should stress that the 
results presented in this paper were obtained in the frame of a more general 
investigation of the interplay of tides and resonance amongst finite mass bodies. 
For this reason, no averaging of the forces was used. Averaged equations are 
fine tools to assess the main effects in restricted there-body problems, but are 
less confident when a third massive body is interacting with the other two. The 
reason is simple: the averaged equations inside and outside resonances are not 
the same (sec Ferraz-Mello, 1987). In a general problem, the semi-major axes 
have secular variations, resonance zones are crossed, one after another, and new 
phenomena are likely to occur which may not appear if the equations are, in 
some way, simplified. 

In view of this, the experiments presented in this work were performed via 
numerical simulations of the exact equations. Although more accurate, this 
approach leads to new problems. Since realistic tidal effects introduce extremely 
slow orbital variations, the integration times are prohibitive. Consequently, we 
have been obliged to enhance these perturbations to be able to obtain significant 
variations in the system, in reasonable integration times. In all our simulations, 
we increased the value of the ratio 

X = 2At{n - n) (1) 

by some 2 or 3 orders of magnitude. Here, At is the tide time lag, H, is the 
rotation angular velocity of the extended body where the tide is raised and 
n the orbital mean motion of the tidally interacting two-body system^. The 
consequence is an acceleration of the evolutionary timescale by the same factor 
increasing x- However, this is not a real scaling of the problem because the 
factor X only enters in one of the right-hand side terms of the equations of 
motion. It is important to keep in mind that the increase of x must not be 
taken too large; the interplay of the migration due to the tide and one resonance 
depends on the speed with which the resonance is crossed, and the probability 
of capture at a crossed resonance lessens as % grows. However, in this paper, we 
are more interested in the phenomena taking place after a capture, which are 
less affected by the enhancement of x as long as this does not introduce large 
forced oscillations about stable stationary solutions. 
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Figure 1: Interaction of an extended body and a mass point. f and g represent 
the resultant of the forces due to tides raised by the mass point on the extended 
body. 



2.1 Tidal forces 

Consider a mass point mi raising tides on an extended body of mass mo. Let 
r be the radius vector of mi relative to mg. The big arrows in figure 1, show 
the central attraction forces between both masses, while f , g are the resultants 
of the forces due to the tidally deformed shape of mg. Since the system of 
forces is in equilibrium they satisfy to the relation f — — g. According with the 
Darwin-Mignard tidal model (Mignard 1981), we have: 

f = -3fc2At ° [2r(r ■ v) + r\rxn + v)] (2) 

where Rq and il are the equatorial radius and the angular rotation velocity 
(vector) of the extended body, k2 is the tidal Love number and v = r. We can 
easily compute the acceleration of the two bodies with respect to an inertial 
frame and subtract them to obtain the relative acceleration: 

„_ G(mo + mi)r ^ (mo + mi)f 



momi 



or 

GMr 



(4) 

r-' m 

where M = niQ + mi is the total mass and m is the reduced mass m = ^mi. 
If we assume that mo 3> mi, then m ~ mi and M ~ mo. 



2.2 The Torque on the Planet 

The momentum of the couple shown in figure is rxf. Since the system is 
isolated, this momentum must be compensated by the momentum of the forces 
acting internally on the extended body with respect to its center of mass: 

T = -rxf. (5) 

is roughly equal to the parameter ~ tane of the Kaula-Goldreich and MacDonald 
(1964) tide models. (Q is the dissipation factor and e is the tide phase lag.) 
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The rotation angular momentum S of the extended body wiU then change fol- 
lowing the law 

S = T (6) 

or, if we assume that the momentum of inertia / of the extended body about 
the rotation axis is constant, 

m = T (7) 

Introducing the expression for the force from equation we obtain 

T = 3fc2A^ ^^ ° [r^n - (r • n)v + vxr] . (8) 

Equations 0-® give the temporal variation of the angular rotation velocity 
of the planet due to the tidal effects. In a series of preliminary simulations 
performed for the Jupiter-Io pair, we found that the perturbation is extremely 
small and the value of remains practically unchanged. As an example, a 
typical run spanning 30 scaled Myrs (corresponding to 1.25 Gyr, considering 
the scaling of x) le&d to a variation of the total orbital angular momentum of 
AS" = +7.4 X 10~^^Mo AU^/yr and, consequently, a variation of the rotation 
period of Jupiter of only +8.8 seconds. 

Following this result, for the rest of our calculations we will consider the 
rotation period of the extended body constant. Of course, a working hypothesis 
of this kind can only be adopted when the relative masses of the orbiting bodies 
to the central one are < 10~^, as in the Jupiter-Io and many planet-exoplanet 
cases. In cases such as the Earth-Moon pair (MacDonald 1964, Touma and 
Wisdom 1994), the mass ratio is much higher 1/80) and the rotation period 
variations should be computed together with those of the orbital parameters. 
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Figure 2: Eccentricities evolution after capture in the 2/1 resonance 
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3 Evolution under Capture 



We now introduce our third body m2 in the system, with an initial orbit exterior 
to TOi. We assume that m2 does not have a tidal interaction with mo and that 
its orbit around too is only perturbed by the gravitational interaction with toi. 
When the semimajor axis of toi increases, mean-motion resonances between 
both satellites are reached. Capture then can take place. We performed a series 
of numerical simulations of the trapping process (and posterior orbital evolution 
within the resonance). The masses and initial conditions used in the simulation 
shown in this section were given by the following parameters: 

mi = 4.684 X IQ-^mo ma = 2.523 x IO^^too 

ai(o) = 5.76i?o 02(0) = 9.40i?o 

ei(o) = 0.0008 62(0) = 0.0003 

Ro = 7.14 X lO^TO n = 1.75 X lO^^s-i 
k2At = 3s 

G = 6.67 X IQ-^^m^kg-^s-^. 

The masses correspond to those of lo and Europa in units of Jupiter's mass. Rq 
and f2 are the equatorial radius and the angular rotation velocity of Jupiter and 
fc2At is a scaled value, approximately 400 times the actual value for Jupiter (the 
estimated Q of Jupiter is very large, see Yoder 1979, Peale 1999). The initial 
semimajor axis of the outer mass, 02(0), corresponds to the current semimajor 
axis of Europa; that of ai(o) is a little less than the current semimajor axis of lo. 
The initial ratio of mean motions (71,1/71,2) was then 2.234. Since the satellites 
were placed above the synchronous orbit (ni < 17), the effect of the tidal friction 
was to increase the orbit of lo. During the whole experiment, ai increases, after 
capture into resonance, from 5.9i?o to 9.8i?o- 

The system evolves with the innermost satellite receding from the planet up 
to the moment where the system is captured into a resonance. In the example 
presented in this section, the 2/1 resonance is soon reached {t ~ 1.2 Myr in the 
scaled time) and the system is trapped in this resonance. Figure 13 shows the 
evolution of the eccentricities after capture. It is worth noting that the variation 
is not nearly monotonic as seen in simplified models. The values increase and 
decrease following the position of the perihelia of the orbits (figure|2|) , a coupling 
that is usually not considered in many studies. The most notorious phenomenon 
is the elbow discontinuities at t ^ 4.65 and t ~ 58 scaled Myrs. The first elbow 
occurs when the two orbits cease having parallel semimajor axes and the second 
elbow corresponds to a return to parallelism after a long evolution through 
stationary solutions with non-aligned apses. In another event, which occurs at 
t ~ 3.05 scaled Myrs, the eccentricity of the inner satellite becomes close to zero 
and the alignment of the semimajor axes change from anti-parallel to parallel. 
(From Am = 180° to Am — 0°.) The several phases of the evolution in low 
eccentricities may be compared to the phase portraits of 2/1-resonant planetary 
systems, given by Callegari et al. (2003). The mode III of Callegari et al has 
a similar behavior showing stationary solutions changing from Am — 180° to 
Am — 0° when a given energy level is crossed. 
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It is worth noting that this behavior is not reported in the Hterature, even 
in complex models such as those studied by Gomes (1998) and Murray et al. 
(2002). We can only assume that their simulations were not extended for a long 
enough time to allow these late stages of the evolution of the system to appear. 
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Figure 3: Top: Evolution of the angles Azu — {w2 — wi), —a\ = Ai — 2A2 + 
VDi and —0-2 = Ai — 2A2 + vj2 after trapping in the 2/1 resonance. Bottom: 
Enlargement of the initial part. 

Comparing those results with the observed configTiration of real bodies, we 
note that the current orbits of lo and Europa have anti-aligned pcriapscs as those 
of the given example before t ~ 3.05. We stress, however, that the time scale 
of our experiment cannot be compared with that of the real satellite system, 
for several reasons. The most important is that our simulations started at a 
time when the bodies were already at the brink of the resonance. Of a lesser 
importance, but significant enough to keep in mind if the real lo-Europa pair is 
meant, is that the adopted time scale flows at least 400 times faster than the real 
one. Additionally, in the actual Galilean system, the Laplacian resonance has 
prevented the lo-Europa subsystem from reaching the deeper libration zone of 
the 2/1 resonance. Thus, it is unlikely that these satellites will ever display such 
turnabouts in libration, even in the far future. For extrasolar planetary systems, 
the current orbits of the two planets orbiting the star Gliese 876 are close to 
have aligned periapses as shown in the above simulations between t ^ 3.05 
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and t ~ 4.65. According to Lee and Peale (2002, 2003), this system may have 
evolved during a certain time in interaction with a remnant dust disk, but the 
evolution stopped when the disk dissipated. As shown by Beauge et al. (2003), 
the actual orbital parameters of this system are somewhat different of those 
corresponding to an exact resonant stationary solution with aligned periapses. 
We may expect that the accumulation of precise observations may give slightly 
different orbital elements or, a second possibility, confirm that the periapses are 
not actually aligned but oscillating about the exact alignment. 

From the dynamical point of view, the more striking behavior is that appear- 
ing in the interval between t ^ 4.65 and t ~ 58 scaled Myrs. During this time, 
the system passes through a sequence of stationary solutions in which the two 
periapses are fixed one with respect to another. However, Aw is no longer 0° 
or 180° as before. In the beginning of this interval. Aw increases very fast from 
0° to about 80°, then, it continues increasing up to 100 — 110° and, thereafter, 
decreases slowly to zero again. These stationary solutions show asymmetric li- 
brations of the angles cri and a2 , a phenomenon until recently only known in the 
restricted asteroidal case (Beauge 1994, Jancart et al. 2002). For illustration 
purposes, figure^ shows, on the right side, planetary orbits in the case of an 
asymmetric stationary solution where the periapses Hi and 112 are separated by 
84°. That figure also shows, on the left side, a symmetric stationary solution 
with the same eccentricities (no matter if not stable). In the symmetric station- 
ary solution with aligned periapses, we have cti = (72 =0, meaning that, not 
only the two periapses are aligned, but the planets have symmetric pericentric 
conjunctions in which both masses pass by the periapses, simultaneously, once 
at each synodic period. The distance of the planets, at the symmetric conjunc- 
tions, is 02(1 — 62) — ai(l — ei), a value that decreases rapidly as 62 grows. In 
the asymmetric stationary solutions, the periapses shift away one from another 
and the conjunctions take place at an intermediate position. This asymmetry 
allows the quantity 02(1 — 62) — ai(l — ei) to approach zero, and even change 
sign, without necessarily leading to an actual collision. Nevertheless, since con- 
junction may occur at places where both planets come very close to each other, 
we cannot expect that the solutions continue to be stable for very large masses. 
Beauge et al. (2003) have shown that these orbits may only exist for planet 
masses less than ^ 10^^ of the central body. 

4 Asymmetric Stationary Solutions 

In the restricted asteroidal problem, it is known that asymmetric solutions only 
exist in the exterior case, that is, when the asteroid is moving in an orbit 
exterior to the planet orbit. Even then, these solutions are only detected in 
resonances of the type p/1, being 2/1 and 3/1 the most important ones (Beauge 
1994). In this section, we search for analogous asymmetric solutions in the case 
of two finite masses, and study their dependence with mi/ 1712. In particular, 
we are interested in detecting the minimum value of the mass ratio for which 
these solutions are still present. According to Hadjidemetriou (2002), symmetric 
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Figure 4: Left: Symmetric stationary solution with aligned periapses. Right: 
Asymmetric stationary solution with |An7| = 84°. The periapses are the points 
Hi and 112. The eccentricities are ei = 0.286 and 62 = 0.3. 

periodic orbits in the planetary 2/1 resonance are only stable for values of 
mi/ 1712 < 1- We will see that asymmetric solutions allow for a larger range of 
masses. 

In these experiments, k2At was increased to 15 s, a value 5 times larger 
than that used in the experiment discussed in previous section. This was a 
limit choice for this parameter; indeed, we can see small fluctuations in some of 
the lines indicating that the variations were not adiabatic enough to keep the 
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Figure 5: Families of stationary solutions of the 2/1 resonance shown in the 
plane (ei, 62). Inner mass is mi = 4.5 x 10^"''. Each curve represents a different 
value of m2; Prom left to right ma = 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.25, 4.5 x 
10~^. The lines corresponding to ma = 4.75, 5.0, 5.25 x 10~^ accumulate at the 
bottom of the set. This accumulation marks the separation between the region of 
asymmetric stationary solutions (above) and symmetrical ones (on and below). 
The diagonal straight line across the figure indicates the limit above which the 
two orbits necessarily intersect. 



9 




90 - 



0,0 0.1 0,2 0,3 0.4 0,5 

ECCENTRICITY e2 

Figure 6: Values of An? for the families of stationary solutions shown in figure|S| 
Top: 1712 > 2.0 X 10~^: the asymmetric solutions bifurcate from Aw = 180°. For 
TO2 > 5.0 X 10^^ no asymmetry appears and Aw does not bifurcate. Bottom: 
TO2 < 1.5 X 10~^: the asymmetric solutions bifurcate from Aw = 180°. 

solutions stationary, forcing small oscillations about the stationary solution. 
4.1 2/1 Resonance 

For this commensurability, the mass mi of the inner orbiting body was fixed at 
4.5 X 10~^ and the outer mass m2 was varied in the range 0.5 — 5.25 x 10^^. 
The results are shown in figures and El 

Figure shows that asymmetric stationary solutions exist for all values of 
mi/m2 above a limit close to 1 (mi/m2 >~ 0.97). As the ratio mi/m2 decreases 
approaching this limit, the interval of eccentricities where asymmetric solutions 
exist decreases to zero. For almost all mass-ratios within this range, asymmetric 
stationary solutions in which the two orbits cross one another exist. The inclined 
straight line in figure El shows the values of ei and 62 such that the apocentric 
distance of mi equals the pericentric distance of m2 . For all initial conditions 
above this curve, the two orbits may intersect. 

The equilibrium values of Aw are shown in figure In the evolution paths 
with m2 > 2.0 x 10^"', the asymmetric solutions bifurcate from symmetric solu- 
tions with aligned periapses; at variance, for m2 < 1.5 x 10^^, the asymmetric 



10 



solutions bifurcate from low-eccentricity symmetric solutions with anti-aligned 
periapses. The two situations are shown in separate figures also because, for 
7712 < 1.5 X 10^^, the eccentricity ei has only a very small variation, while the 
62 has a large variation and thus the variations are more clear shown if it lies 
in the horizontal axis. 

It is worth comparing these results on the 2/1 resonance with the periodic 
orbits of Hadjidemetriou (2002). His conclusion that stable solutions are only 
found with mi jmi < 1 is true only if we restrict the domain to large eccentricity 
symmetric orbits. In the general case, stable stationary solutions can be found 
for much higher values of the mass ratio. 

One last point to mention is that the asymmetric solutions appear as a 
bifurcation point where symmetric solutions change from stable to unstable. 
This bifurcation generates two distinct families of asymmetric solutions, each 
independent of the other. When a trapping occurs, the system may be captured 
in either one or the other of the two new centers. This is clearly seen in figures 
where approximately half of the solutions departed in one direction and half 
in the opposite one. 

4.2 3/1 Resonance 

In the restricted asteroidal problem, asymmetric solutions exist for all exterior 
resonances of the type p/1 (i.e. 2/1, 3/1, 4/1, etc.). A series of simulations 
was then done to see whether the general three-body problem also presented 
asymmetric points in the 3/1 resonance. As in the previous subsection the 
relative mass of the innermost orbiting body was fixed at 4.5 x 10~^ and the 
outer value was varied from 2.5 x 10^'°^ to 5.5 x 10~^ in steps of 0.5 x 10~^. The 
remaining parameters were chosen as: 



The results are shown in figs. [7| and |S1 The capture into resonance occurs 
quickly and Aro reaches 180° remaining there up to the bifurcation and switch- 
ing to an asymmetric stationary solution. Periodic orbits with aligned periapses 
were not seen, although, in the very beginning of the simulations Atu appeared 
temporarily oscillating about 0° with a large amplitude. The more characteris- 
tic feature of the 3/1-resonant asymmetric stationary solutions is the boundary 
between symmetric and asymmetric regions at 62 0.11 (elbows line in figure 
Uj). In this resonance, asymmetric stationary solutions exist for all values of the 
mass ratio (greater or smaller than unity) up to the limit corresponding to a 
curve whose elbows occurs for 6i ~ 1. For mass ratios larger than this, the 
stationary solution lines will remain below the 62 ^ 0.11 limit for all 61 < 1. 



ai(o) = 5.85i?o 
6i(o) = 0.0010 



02(0) = 12.4i?o 
62(0) = 0.0004 
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Figure 7: Families of stationary solutions of the 3/1 resonance shown 
in the plane ei,e2; mi = 4.5 x 10~^; from left to right m2 = 
0.5, 1.0, 1, 5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5 x 10"^ The line formed by the el- 
bows of the curves, at 62 ~ 0.11, separates the symmetric (below) and asym- 
metric (above) stationary solutions. 
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Figure 8: Values of Am for the families of stationary solutions shown in figure 
dfor which m2 > 2.5 x 10^'^; In the 3/1 resonance the asymmetric solutions 
bifurcate from symmetric solutions with anti-aligned periapses (Atn = 180°) 

5 Conclusions 

In this paper, we presented a series of numerical simulations of the evolution 
of systems of two massive orbiting bodies after capture into the 2/1 and 3/1 
resonances. Albeit with a change in the timescale of the solutions, the same 
behavior can represent the evolution of a pair of satellites around one planet 
(the data used were those from the lo-Europa pair) or two exoplanets orbiting 
close to a star. The simulations were done using tides in the central body as 
the source of the non-conservative perturbation, but the results do not depend 
on the particular force used, and should happen in similar way with other non- 
conservative forces provided its effects are sufficiently slow to have adiabatic 
variations. Thus, the same behavior is expected to hold in the evolution of exo- 
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planets captured into a 2/1 or 3/1 resonance due to (for example) interactions 
with a residual disk of matter or any other non-conservative force. For instance, 
Lee and Peale (2003) have obtained results very similar to those discussed in 
section 3 in a simulation of a system of 2 planets in which the outer planet is 
driven inward by torques exerted on it by outside nebular material. 

Finally, we have shown that the evolution of resonant trapped massive bod- 
ies is not simple. Our results show a surprising richness of solutions: symmetric 
librations around Aoj = and around Aoj = tt, asymmetric solutions with 
a wide range of values of Aro, turnabouts from one configuration to another 
during the secular variation of scmimajor axes, etc. Many of these solutions 
exist for large values of the eccentricities. As our knowledge of extrasolar plan- 
ets continues to grow, it will be interesting to see how many of these possible 
configurations actually occur in the real world. 
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